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Abstract. The interpretation of water lines in red supergiant stellar at- 
mospheres has been much debated over the past decade. The introduc- 
tion of the so-called MOLspheres to account for near-infrared "extra" 
I I absorption has been controversial. We propose that non-LTE effects 

Ph should be taken into account before considering any extra-photospheric 

^0 contribution. 

^ After a brief introduction on the radiative transfer treatment and 

f^ , the inadequacy of classical treatments in the case of large-scale systems 

I such as molecules, we present a new code, based on preconditioned 

O Krylov subspace methods. Preliminary results suggest that NLTE ef- 

■^^ fects lead to deeper water bands, as well as extra cooling. 

1 Water in RSG atmospheres 

^^ The interpretation of water lines in red supergiant stellar atmospheres has been de- 



^H bated for more than a decade. For example, Tsuji (2000) shows that near-infrared 
^-H water bands observed with the Infrared Satellite Observatory (ISO) low resolution 
"^ spectrograph toward fi Cep appear significantly deeper than expected from a pure 
"^ photospheric, ID, local thermodynamic equilibrium (LTE) and hydrostatic syn- 
^^ thetic spectrum. In order to quantify this extra absorption, he considered the pos- 
sibility to add an ad hoc quasi-static molecular shell baptized " MOLsphere" . But 
what was initially introduced as toy model gained support through interferometric 
observations which reveal a larger spatial extension in water bands compared to 



m 



the stellar radius determined in the continuum (see Perrin et at. 2005 Ohnaka 
[20041 ). 



^ However the MOLsphere model is still subject to some weaknesses. First of all, 

no model is up to now able to explain the formation of such extra-photospheric 
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layer (even gas levitation by shock wave cannot account for the parameters such 
as density and water abundance proposed in MOLsphere models). One should 
also keep in mind that, due to the rather poor coverage of the Fourier plane by 
current infrared interferometers, the interpretation of such observations strongly 
rely on a presupposed model. More importantly, based on high resolution spectra 
obtained with TEXES,|Ryde et ar|(|2 006[ ) showed that water lines around 12 /im, 
predicted in emission by MOLsphere models, are in fact observed in absorption 
and are probably of photospheric origin. 

In order to reconcile the interpretation of observations in various spectral 
ranges, we want to test how important non LTE (NLTE) effects are and whether 
they may mimic a MOLsphere through a strengthening of absorption lines. These 
NLTE effects may also affect the thermodynamical structure of the atmospheres 
through e.g. the contribution of molecules to the cooling function and the ra- 
diative pressure. If it leads to a more extended atmosphere and a shift of the 
water line contribution function, this could potentially explain the interferometric 
observations. 

2 General context of radiative transfer in molecular lines 

2.1 Basics of radiative transfer 

In order to determine the specific intensity /^ in a line, one has to solve the 
radiative transfer equation (RTE). For a ID plane parallel geometry, this equation 
is a simple ordinary differential equation of the first order, 

fi^=L-S,. (2.1) 

Solving the RTE under NLTE conditions is however not straightforward, because 
of an intricate non-linear and non-local coupling between level populations and 
radiation. Indeed the determination of the source function for a given line due to 
a transition i ^ j at an optical depth r requires the knowledge of the populations 
on the levels i and j at this optical depth. 

^^.M = ^i-l' (2-2) 

where Xk = Tik/ Qki with Uk being the relative population of level k and gk is its 
statistical weight. The other symbols have their usual meanings. These popu- 
lations can be derived by solving the statistical equilibrium system of equations, 
which is a set of rate equations for each energy level. 

N N 

Yl ^^i + Yl i^ij^ij + ^coiCij (T)) 

j<i j^i 

N N 

^ UjAji + ^ Uj {BjiJji + ricoiCji (T)) | = 0, 



(2.3) 
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where J ji denotes the mean radiation field for the i ^ j transition, integrated over 
the line profile. This system has in general no analytical solution. However, if one 
considers that the probability of radiative (de) excitation is very small compared 
to the collisional (de)excitation, either because the radiative field (through J ji) is 
weak or because the density of colliders ricoi is large, the system can be simplified 
and leads to an analytical solution for the level populations, being the Boltzmann 
distribution. The populations then only depend on the local kinetic temperature 
at r. This is the LTE hypothesis. The source function is simplified to the Planck 
function which only depends on the local temperature as well. On the contrary, if 
the radiative field is strong and/or the collider density is low, one has to consider 
the radiative terms associated to J and invert this system of equations. The main 
difficulty is that the mean intensity field J depends on the solution of the RTE: 

Jij{r) = - Mr) L{r,^)d^du (2.4) 

This explicit dependence on the solution means that to solve the statistical equi- 
librium in order to determine the source function, one need to know the intensity 
at different locations and for different frequencies simultaneously. The radiative 
transfer is thus a non-local problem where the physical state at a given r is coupled 
with every other layer. Furthermore, for a multi-level system, the problem is also 
non-linear. Because of this coupling, the source function and the mean intensity 
field are usually determined iteratively. This is the classical way to solve NLTE 
radiative transfer. 

2.2 The water molecule 

Solving the statistical equilibrium system requires molecular data: the energy and 
statistical weight of each level, the Einstein coefficients Aij (including the selection 
rules), and finally, the collisional rates Cij. For the water molecule, we have used 
the line list and radiative coefficients from Barber et aL\ (|2QQ6). The radiative 



transitions we have taken into account are displayed in Fig. M The collisional 



rates were computed by Faure & Josselin (2008) thanks to an original method 



of extrapolation for high J (relevant to the temperatures met in red supergiant - 
RSG - atmospheres). 

A first rough estimate of the deviation from LTE can be obtained by evaluating 
the competition between radiative and collisional processes. A good indicator of 
this deviation is thus the ratio between radiative and collisional de-excitation rates, 
that is, the ratio between the critical density and the density of colliders for each 
energy level. 

,(,) <"(^) 1 ^^^^^ (2 5) 



n^ 



\t) n-ol{T)Y!-}^Cij{T{T)) 



If this ratio rji is much smaller than 1 for every level i, collisions dominate ev- 
erywhere over radiative transitions and the solution of the statistical equilibrium 
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Fig. 1. Grotrian diagram of the water molecule with the radiative transitions taken into 
account in this work. Jn is the total angular momentum. 



equations tends toward the Boltzmann distribution (LTE). On the contrary, if rji 
is greater than one, NLTE conditions are met. 
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Fig. 2. Ratio between the critical density and the density of colliders (e and H2) in a 
Betelgeuse-like atmosphere. The solid line indicates the transition rfi — 1. 

This ratio is computed and displayed in Fig. [2] in the case of water in a typical 
RSG atmosphere, rji is much larger than 1 in the upper atmosphere, i.e. where 
the water line contribution function is expected to peak. This results mostly from 
the water abundance profile and the excitation temperature. An accurate NLTE 
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treatment of the RTE is thus required to properly analyze water spectra. 

2.3 The weakness of classical methods for large-scale systems 

The classical methods generally used to solve the RTE rely on the operator split- 



ting technique, which is derived from the Lambda Iteration (LI) (Trujillo Bueno 
fc Fabiani Bendicho||1995). The LI method consists in iteratively computing the 



Schwarzschild-Milne integral with the use of an integration operator A. A inte- 
grates the source function over the domain (i.e. the stellar atmosphere) in order 
to determine the mean intensity field. However, LI is inefficient and inaccurate 
in many cases, as it does not converge toward the true solution. This is due to 
the fact that the spectral radius of the amplification matrix of this scheme A^^^ 
depends on the NLTE parameter e (i.e. the photon destruction probability; the 
lower e, the stronger the NLTE effects) and the total optical depth of the medium 
T, as A^^^ c^{l-e){l-T-^)-^. It is close to one if e is smah (i.e. strongly NLTE) 
and/or if T is large. In order to improve its convergence toward the true solution, 
a better conditioning can be obtained through operator splitting methods, such 



as ALI (Jacobi), Gauss-Seidel or Successive OverRelaxation method SOR (Tru- 
jillo Bueno fc Fabiani Bendicho||1995), which reduce the spectral radius (see Fig 



3|. Unfortunately, for large discretized and strongly NLTE problems, the spectral 
radius remains large, and the convergence slow (e.g. typically ~ 100 iterations 
in case B in Fig. [3|. This slow convergence becomes problematic if the model 
atom or molecule is large, as for each iteration the large statistical equilibrium 
system has to be solved. For water, this system is huge (about 1000 energy levels, 
corresponding to about 15,000 radiative rates and more than 330,000 collisional 
rates; see Fig. fl]) and the total time needed to converge becomes prohibitive. Fur- 
thermore, the accuracy of these methods may be poor (see Chevallier et a/.|[2003 



especially Fig. 9). This is intrinsically due to the static nature of the iterations. 

3 New method and preliminary results 

3.1 Methodology 

Rather than determining iteratively the source function, the new method described 
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(2012) is to determine, in a unique step, the RTE and the populations 
in every layer of a discretised model atmosphere. This is done by inverting a 
generalised multi-zone statistical equilibrium system with a non-linear solver based 
on non-static methods such as sub-space Krylov methods ( |Saad||2003| . 

The first step consists in rewriting the RTE with an explicit dependence on 
populations which become the new unknowns. Moreover, because of the non- 
locality of the RTE and the need to couple the sets of statistical equilibrium in 
every layer, we use a unique transition-independent scale (e.g. ry = r^oo)- 



dry 



fi^ = au{xj,Xi)Iu -ju{xi) (3.1) 
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Fig. 3. Spectral radius (left) and convergence properties (right) for two static iterative 
methods (LI: upper part; ALI with A"^ = diagonal of A: lower part). For each method, 
the convergence is shown for two spectral radii, marked with letters A (upper panel) and 
B (lower panel) in the right-hand panels (the true solution is given there by the orange 
line) . 



Then, the formal solution of this equation can be determined by the use of its 
Green function and its integration over the model atmosphere. By including this 



solution into Eq. 2.4, an analytical form of the mean intensity is obtained, with 
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an explicit dependence of the population: 



J ij 



_out 

(rv) = / fij{xi{t)) ' Kij[xi,Xj]{t,TY)dt (3.2) 

Jrkf 
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where fij {xi{t)) denotes the explicit dependence on populations (see 

|2Q12[ for more details) and Kij is the kernel of the RTE. This expression is included 



into the statistical equilibrium equations of each layer (labelled by /c, 1 < /c < P) 
and is rewritten in a matrix form. 



M^ 



J^(x^x^...,x^)l .x^ = (3.3) 



-k 

Because the RTE is non-local, J requires an integral over ry, and thus the knowl- 
edge of the populations everywhere (x^,x^, ...,x^). The well-adapted unknown 
vector is thus X = (x^,x^, ...,x^) and the full closed system is built: 

/ Ml (X) \ / xi \ / \ 

r(x).x= •. • i =h (3-4) 

\ M^ (X) / \ x^ / \ / 

Solving the RTE and the statistical equilibrium equations then consists in finding 
the root of: 

r(X) • X = F(X) = (3.5) 

The inversion of this system gives the populations everywhere in the model which 
enables the computation of the source function. The emergent spectrum and any 
other physical quantity are easily post-processed. An interesting advantage of this 
approach is the analytical form of this N x P integral equation (for N energy levels 
and P atmosphere layers), and thus the possibility to differentiate each equation 
with respect to populations on every level and every layer in order to build an 
analytical Jacobian matrix (which is sparse) very rapidly. 

3.2 A new MPI code: MOrad 

This method has been implemented in a fully parallelised MPI code which takes 
as input a model atmosphere and molecular data. The chosen MPI strategy is a 
domain decomposition on spatial cells. The code computes the multi-D function 
F(X) and its Jacobian Matrix. It is interfaced with the efficient non-linear and 
scalable library PETSc, which solves the equation with the "generalized minimal 
residual method" (GMRES) method, well adapted for large and sparse problems 
after preconditioning of the system (]Saad fc Schultz||l986). 



3.3 Preliminary results 

The departure coefficients in a Betelgeuse-like model atmosphere are displayed 
in Fig. [4j As expected from the critical densities, the strongest deviations from 
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LTE are met in the upper atmosphere. Furthermore, the deviation from LTE is 
all larger when the energy is higher. The dispersion of the departure coefficients 
within a given vibrational level is remarkably low for highly-excited vibrational 
bands. In other words, these bands are under-populated with respect to LTE, but 
rotational levels are in relative LTE within these bands. This partially validates 
the potential use of the super-level technique for water (see Schweitzer et a/.|2QQQ ). 
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Fig. 4. Departure coefficients in a RSG model atmosphere. Vibrational bands are coded 
with different colours. 

As shown in Fig. |5]the emerging spectrum exhibits water bands which are re- 
markably deeper compared to LTE predictions. NLTE seems thus to qualitatively 
mimic a MOLsphere. 



4 Conclusions and perspectives 

Our preliminary results show an important deviation from LTE in water lines, 
which should be taken into account in the interpretation of infrared spectra of 
RSG and especially when discussing the MOLsphere hypothesis. These NLTE 
effects may also have strong implications on the thermodynamic structure of RSG 
atmospheres. For example, preliminary calculations show an increase by up to a 
factor of 4 in the cooling due to water molecules. The structure and the extent of 
these atmospheres may thus be strongly affected. For example, a shift outward of 
the contribution function peak is possible, which would impact the interpretation 
of interferometric observations. Finally, radiative pressure on molecules may play 
a role in the mass loss process (see Josselin fc Plez||2QQ7). 
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Fig. 5. Sample NLTE (red) and LTE (pink) water spectra emerging from a RSG. 
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